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We  consider  positivity  preserving  property  of  first  and  higher  order  finite  volume  schemes 
for  one  and  two  dimensional  compressible  Euler  equations  of  gas  dynamics.  A  general  frame¬ 
work  is  established  which  shows  the  positivity  of  density  and  pressure  whenever  the  under¬ 
lying  one  dimensional  first  order  building  block  based  on  an  exact  or  approximate  Rie- 
mann  solver  and  the  reconstruction  are  both  positivity  preserving.  Appropriate  limitation 
to  achieve  high  order  positivity  preserving  reconstruction  is  described. 
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1  Introduction 


To  solve  a  scalar  conservation  law 


Ut  +  div  /(«)  =  0  (1.1) 

with  possibly  discontinuous  solutions,  one  usually  studies  the  total  variation  stability  (in  one 
space  dimension)  or  L°°  stability  (in  multi  space  dimensions)  of  the  numerical  schemes.  The 
schemes  can  usually  be  generalized  to  systems  via  local  characteristic  decompositions  and 
usually  work  equally  well  numerically.  However,  no  stability  property  can  be  automatically 
carried  over  to  the  nonlinear  system  case.  For  example,  most  second  order  TVD  schemes,  or 
even  some  first  order  monotone  schemes,  when  generalized  to  compressible  Euler  equations 
of  gas  dynamics,  do  not  always  preserve  the  positivity  of  density  and  pressure.  This  may 
cause  problems  in  practical  calculations  when  the  solution  is  near  vacuum,  for  example  in 
the  computation  of  blast  waves. 

When  considering  the  compressible  Euler  equations  of  gas  dynamics,  a  natural  stability 
condition  for  the  numerical  approximation  is  positivity  preserving  for  density  and  pressure. 
Another  possible  a  priori  bound  is  the  maximum  principle  for  the  specific  entropy  (Tad- 
mor  [10]),  which  seems  extremely  difficult  to  preserve  for  second  and  higher  order  schemes 
(Khobalatte  and  Perthame  [5]). 

The  positivity  of  density  and  pressure  is  already  interesting  because  it  allows  one  to  derive 
a  rigorous  CFL  condition  (limitation  of  the  time  step  to  use).  Since  we  consider  conservative 
schemes,  it  also  allows  one  to  obtain  a  priori  L1  bounds  on  the  density,  momentum  and  energy. 
In  this  sense,  positivity  of  density  and  pressure  is  a  weak  stability  condition.  Of  course  for 
these  nonlinear  systems,  such  a  weak  stability  is  not  enough  to  assert  the  convergence  of  the 
method  and  estimates  on  the  derivatives  are  usually  needed  for  that  purpose.  Note  however 
that  for  a  large  class  of  2  x  2  systems  in  one  dimension,  weak  stability  (say  in  L°°)  is  enough 
to  prove  the  convergence  of  the  method. 

In  this  paper  we  will  provide  a  general  framework  and  illustrate  by  several  examples  the 
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way  to  impose  positivity  of  density  and  pressure  for  finite  volume  schemes,  for  one  and  two 
space  dimensions  and  for  first  and  higher  order  accuracy.  Of  course  the  first  ingredient  is  a 
positivity  preserving  exact  or  approximate  Riemann  solver,  such  as  Godunov,  Lax-Friedrichs, 
Boltzmann  type  (Perthame  [6]),  and  Harten-Lax-van  Leer  [3].  For  a  first  order  scheme,  this 
is  enough  to  guarantee  positivity  of  density  and  pressure,  even  for  two  dimensional  problems 
set  on  an  arbitrary  triangulation.  However,  for  higher  order  finite  volume  schemes,  the 
reconstruction  must  also  satisfy  such  a  property.  We  will  show  that  only  the  nodal  values 
needed  for  constructing  the  flux  along  a  cell  edge  must  be  positive,  in  order  to  obtain  a 
positive  scheme.  The  reconstructed  function,  which  is  a  piecewise  polynomial  and  can  be 
obtained  in  ENO  spirit,  needs  not  be  positive  everywhere. 

This  paper  is  organized  as  follows:  we  first  prove  a  positivity  result  and  apply  it  to  a 
third  order  scheme  in  one  space  dimension  in  Section  2.  Then,  in  Section  3,  we  consider 
first  order  positive  schemes  in  two  space  dimensions  for  arbitrary  triangulations.  The  fourth 
section  is  devoted  to  second  order  schemes  in  two  space  dimensions.  In  the  Appendix  we 
recall  why  Lax-Friedrichs  scheme  is  positivity  preserving. 

2  Positivity  in  One  Space  Dimension  and  a  Positive 
Third  Order  Scheme 

In  this  section  we  consider  the  one  dimensional  compressible  Euler  equations  for  perfect  gas: 

Ut  +  F{U)X  =  0,  t  >  0,  x  €  R  (2.1) 

with 

U  =  (p,pu,E),  E=)^pu2  +  pe  (2.2) 

F{U)  =  (pu,pu2  +p,(E  +  p)uj  ,  p=(~/-l)pe  (2.3) 

where  p  is  the  density,  u  is  the  velocity,  E  is  the  total  energy,  p  is  the  pressure,  e  is  the 
internal  energy,  and  7  >  1  is  a  constant  (7  =  1.4  for  air). 
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Let  us  recall  for  later  purpose  that  for  the  pressure  law  in  (2.3),  the  speed  of  sound  is 
given  by  c  —  ^7(7  —  1  )e  and  thus  the  three  eigenvalues  of  the  system  are  u,u  ±  c. 

We  give  conditions  on  numerical  schemes  in  order  to  satisfy  the  positivity  property,  by 
which  we  mean  that  the  resulting  value  of  U  should  satisfy  p  >  0(P]  >  0  on  the  discrete 
level)  and  p  >  0  (p”  >  0  on  the  discrete  level),  if  the  initial  condition  satisfies  those  positivity 
conditions. 

We  first  present  the  finite  volume  schemes  under  consideration,  then  we  give  a  general 
positivity  theorem,  and  finally  we  show  how  a  third  order  reconstruction  can  be  modified  in 
order  to  satisfy  the  assumption  of  the  general  theorem. 

A  general  finite  volume  scheme  can  be  written  as 

U*+l  =  (7"  -  A  [/,  (t/;+i,(/-+})  -  h  {uik,u-_ j)]  (2.4) 

where  n  >  0  refers  to  the  time  step  At  and  j  €  Z  to  the  uniform  space  discretization  of 
the  size  Ax  (to  simplify  the  presentation),  and  A  =  U™  are  approximations  to  the  cell 
averages  of  (I  in  the  cell  Cj  =  (a^.i, a^+i)  at  time  level  n,  and  U~+i,  U++ 1  are  high  order 
approximations  of  the  nodal  values  Un(xj+i)  within  the  cells  C,  and  C,-+ 1,  respectively.  These 
values  can  either  be  evolved  as  independent  degrees  of  freedom  such  as  in  the  discontinuous 
Galerkin  finite  element  method  (see,  e.g.,  [1]),  or  be  reconstructed  from  the  cell  averages 
If-.  Let  us  recall  that  the  general  ENO  philosophy  [4]  allows  us  to  reconstruct,  from  the  cell 
averages  U ]  ,  a  piecewise  polynomial  function  Un(x)  which  is  high  order  accurate  (r-th  order 
if  Un(x)  is  piecewise  (r  —  1  )-th  order  polynomials),  and  conserves  the  local  mean: 

h  lc,u'{x)ix  =  Tr>  (2-5) 
The  nodal  values  needed  in  scheme  (2.4)  can  then  be  set  to 


since  the  function  Un(x)  is  discontinuous  at  the  cell  interface  Xj+ 1. 
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What  remains  to  be  explained  for  scheme  (2.4)  is  the  flux  function  h(U,  V).  This  is 
assumed  to  be  an  exact  or,  in  most  cases,  an  approximate  Riemann  solver.  In  particular, 
h(U,  K)  is  a  Lipschitz  continuous  function  of  both  arguments,  and  is  consistent  with  the 
physical  flux  F(U)  in  (2.3):  h(U,  U)  =  F( U).  We  will  further  make  the  following 

Assumption  1:  h(-,  •)  produces  a  one  dimensional  first  order  scheme  which  satisfies  the 
positivity  property  under  a  CFL  condition: 

A||(|«|  +  c)|L<a0  (2.7) 

i.e.,  if  p"  >  0  and  p"  >  0  for  all  j  €  Z ,  then  the  solution  £/"+1  of 

(/»+'  =  (/»  —  A  [h(U«„ ,  U?)  -  h(U",  CJL, )]  (2.8) 

also  satisfies  p"+1  >  0  and  p”+1  >  0  under  the  CFL  condition  (2.7). 

We  will  call  such  a  h(-,  -)  positivity  preserving  under  the  CFL  condition  (2.7).  Examples 
of  positivity  preserving  (approximate)  Riemann  solvers  include  Godunov,  Lax-Friedrichs  (see 
the  Appendix),  Boltzmann  type  [6]  and  Harten-Lax-van  Leer  [3]. 

Our  first  result  is 

Theorem  1.  Assume  A(-,-)  satisfies  Assumption  1.  If  the  reconstructed  nodal  values  U.  t 

j+2 

have  positive  density  and  pressure  for  all  j  6  Z,  then  the  full  scheme  (2.4)  is  positivity 
preserving  under  the  CFL  condition 

AIKIul  +  c)!^  <  oa0  (2.9) 

where  0  <  a  <  1  is  sufficiently  small  such  that 

Pj  -  a(uJ+l  +  ut-\)  '■=  V>  (2-10) 

has  positive  density  and  pressure  for  all  j  €  Z. 

Remark  1.  Of  course  a  =  0  works.  However,  since  U±  L  are  close  to  TF-,  we  can  expect 
that  a  is  close  to  j. 
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Remark  2.  Here  we  restrict  ourselves  to  first  order  Euler  forward  in  time.  TVD  type  high 
order  Runge-Kutta  time  discretization  (Shu  and  Osher  [9])  will  keep  the  validity  of  Theorem 
1,  as  it  was  shown  in  [5]. 


Proof  of  Theorem  1:  Let  us  introduce  three  “first  order”  schemes  within  the  cell  Cy. 


Vj 

VJ 


v3-\ 

vu- 

J  2 


sK^.<H(V')] 


(2.11) 


The  three  of  them  are  of  the  type  (2.8),  with  possibly  £  in  the  place  of  A.  Therefore  under 
the  CFL  condition  (2.9),  the  values  f/+,  U~  and  V  all  have  positive  density  and  pressure. 
We  set 

(7”+'  =a(Uf  +  £r)  +  V,  (2.12) 

Then  the  linear  combination  of  the  above  three  equalities  in  (2.11)  with  weights  (o,  l,a)  just 
gives  the  scheme  (2.4),  thanks  to  the  definition  of  Vj  in  (2.10). 

Finally,  (2.12)  implies  that,  by  concavity,  U”+1  has  also  positive  density  and  pressure. 
The  pressure  is  indeed  a  concave  function  of  ( p,pu ,  E). 


□ 


Theorem  1  tells  us  that,  if  we  use  a  positivity  preserving  approximate  Riemann  solver, 
then  we  only  need  to  worry  about  positive  density  and  pressure  in  the  nodal  values  (/*  , , 
either  from  the  reconstruction  or  from  direct  evolution.  This  can  usually  be  achieved  by  a 
further  limitation  upon  U±  L,  such  that  positivity  of  density  and  pressure  is  enforced  and 
accuracy  is  preserved  in  smooth  regions.  We  now  show,  as  an  example,  how  this  can  be 
done  for  a  third  order  reconstruction.  Let  us  consider  a  typical  cell  (—6,6)  where  6  =  —,  in 
which  we  are  given  three  cell  averages  (p,pu,E)  with  positive  density  and  pressure.  Let  us 


build  three  quadratic  polynomials 


x * 

p(x)  =  Po +  />i*  + Pay 

x2 

u(x)  =  u0  +  «ix  +  u2—  (2.13) 

/  \  x 2 
p(x)  =  Po+pix-fp2y 

2  — 

obtained  by  (i)  fitting  the  average  of  (p,pu,  +  ^y)  to  {p,~pu,  E )  and  (ii)  fitting  the  nodal 

2 

values  (p,pu,  +  ^£y)(±6)  with  the  values  obtained  by  a  third  order  reconstruction,  say 
from  the  ENO  procedure.  These  might  be  associated  to  a  negative  density  and  pressure. 
Then  we  need  an  additional  limitation.  It  can  always  be  performed  easily  as  shown  in 

Proposition  2:  Given  the  three  quadratic  polynomials  in  (2.13),  we  can  always  perform 
a  limitation  of  their  coefficients  so  that  (i)  the  means  of  the  three  conserved  variables  are 
preserved,  and  (ii)  the  nodal  values  have  positive  density  and  pressure. 

Remark:  The  precise  limitation,  and  the  deduction  of  the  appropriate  coefficients  in  (2.13) 
are  given  in  the  proof  below. 

Of  course,  combining  Theorem  1  and  Proposition  2,  we  obtain  a  positive  preserving 
scheme.  A  similar  construction  for  second  order  schemes  can  be  found  in  [6]  and  [5]. 

Proof  of  Proposition  2:  First  of  all,  let  us  fix  some  notations.  Being  given  the  cell  averages 
(p,~pu,  E),  we  also  define  u  and  p  through 

pu  =  pu,  E  =  (2.14) 

2  7  —  1 

Notice  that  we  are  given  p  >  0  and  p  >  0. 

x2 

First  step:  The  limitation  on  the  density  is  rather  easy.  Since  p  =  P0+P2 y,  we  can  modify 
pi  without  altering  the  conservation  of  mass.  To  enforce  positivity  of  p(±£)  is  equivalent  to 
impose  a  simple  limitation  on  p\i 


But  the  right-hand-side  of  (2.15)  could  be  negative.  Therefore  we  first  impose,  for  instance, 


5p  .  3  p 

Y  Po  T' 


P  —  PO  +  p2 ' 


(2.16) 


Indeed,  this  implies  that 


|p2|y  =  3|p2|y  =  3| p  -  p0\  <  ~  <  po 


and  now  (2.15)  can  be  imposed  too. 

Second  step:  Next  the  conservation  of  momentum  implies  the  following  relation  which 
gives  uQ  in  terms  of  ut,  u2  and  p(x): 

82  8 2  8 4 

pu0  =  pu  -  P\u\—  -  u2{p0  —  +  P2^)  (2.17) 

And  the  conservation  of  energy  gives,  after  some  easy  calculations 

_  ,  1  ,  ,2  ,  (8*  84\  8 4  286  82 

P  =  { 7-1)  p(u°~u)  +  u\\poj+P2~j+uiu2P\j-  +  p2U2  —  +P0  +  P2—  (2.18) 

As  in  the  density  case,  we  will  impose,  in  order  to  get  a  positive  pressure  at  the  nodes,  a 
simple  limitation  on  pj 


|Pl|$  <  PO  +P2- 


(2.19) 


and  thus  we  need  some  limitations  in  order  to  guarantee  that  the  right-hand-side  of  (2.19) 


is  positive.  As  before  we  first  impose  the  limitations  on  po,ui,u2: 


Po<-P 


«i(7“l)  Pov  +  IPif 


“2(7  -  I)  fPii 


.  <$3  84\  p 

'  l/?llTo  +/?2!oJ  <  I2 

.85  86\  p 

*10  +  56  /  <  12 


(2.20) 

(2.21) 


(2.22) 


This  defines  a  unique  uq  through  (2.17)  and  a  unique  p2  through  (2.18).  For  these  values  we 

deduce  from  (2.18),  using  (2.20): 

P  ^  _  2 

6  < 

,  1.,  _,2  ,  2  /  82  ,  84\  84  286  1  (  82\ 
=  (7-1)  ~2 P{Uo-u )  +u1lpoj  +  P2  —  j+pxU1U2—  +  p2u'2—  +  -IP0+P2-7H 

/  /  ,  J  2  (  P  ,  ,<53\  2  /  86  .  ,8*W  1  (  t2\ 

<  (7-1)  [u,  («, J  +  n-  +  J  +  (p2-  +  l/»l]o)J  +  3  (Po  n  j) 
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and  from  (2.21)-(2.22)  we  deduce  that  p0  +  p-i  T  >°  and  thus  (2.19)  can  also  be  imposed. 

This  concludes  the  description  of  a  possible  set  of  limitations,  which  proves  Proposition  2. 
Notice  that  these  limitations  (2.15),  (2.16),  (2.19),  (2.20),  (2.21)  and  (2.22)  just  amount  to 
avoid  gradients  or  second  order  derivatives  of  the  order  ^  (or  larger)  and  thus  they  preserve 
the  accuracy  in  smooth  regions. 


3  First  Order  Positive  Schemes  in  Two  Dimensions 


We  now  consider  the  equations 


Ut  +  div  F(U )  =  0  t  >  0,  x  €  R2 


where 


U  =  (p,pu,E),  ueR2,  E  = -p\u\2  +  pe  (3.2) 

F(U)  =  (pu,pu®u  +  pI,(E-\-p)u),  p  =  (7  —  l)pe  (3.3) 

We  consider  finite  volume  schemes  set  on  a  triangulation  C.  The  control  volumes  will  be 
the  triangles  K  £  C.  For  each  triangle  K  we  denote  by  (/£) i<q<3  the  length  of  its  three  edges 
(e£-)i<a<3,  with  outward  unit  normal  (^)i<a<3-  Finally,  K(at)  will  denote  the  neighboring 
triangle  along  and  |A'|  the  area  of  the  triangle  K.  Then,  we  consider  the  scheme 

A  t  ** 

u’k"  =  U’K  -  £  k(i/Sw,i/*,  v%)rK  (3.4) 

for  some  (approximate)  Riemann  solver  h(U,  V ,  v)  in  the  direction  v.  We  recall  some  classical 
examples  in  the  Appendix.  The  basic  properties  of  h  are  now 

h{U,V,v)  =  —h(V,U,—i/),  ( conservativity )  (3.5) 

h(U,U,v)  =  F{U)  ■  v  ( consistency )  (3.6) 

Next,  we  also  impose  a  positivity  condition  for  the  one  dimensional  solver  obtained  by 
fixing  v  (see  the  Appendix  for  examples).  This  is  the  same  as  Assumption  1  in  the  previous 
section,  which  now  means  that  the  solution  (/,  a  four  component  vector,  of 


U  =  U  -  A  [h(V,  U ,  u)  -  h(U ,  W,  v)} 
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has  positive  density  and  pressure  as  soon  a s  U,V,W  do.  Now,  h  is  a  four  component  vector. 
The  main  result  of  this  section  is 

Theorem  3:  Let  h(-,  •,  •)  satisfy  (3.5)-(3.6)  and  the  one  dimensional  positivity  property  As¬ 
sumption  1.  Then  the  scheme  (3.4)  satisfies  the  positivity  property  under  the  CFL  condition 

/£A£(|uk|  +  cK)  <  c*o|A'|,  for  all  K  £  C  (3.8) 

Of=l 


Proof  of  Theorem  3:  We  define 


Ua  = 


la 


At 

S3.,  \I°K\  K  “  w 


<K  [M (3.9) 


Since 


E  iWU"K,u"K,*i)  =  f(U"k)  ■  E  =  o 


a=l 


Qf  =  l 


we  have  U^+l  —  5Za=i  Ua •  we  conclude  as  in  the  one  dimensional  case. 


4  Second  Order  Positive  Schemes  in  Two  Dimensions 

We  extend  the  result  of  the  previous  section  to  second  order  schemes,  i.e.,  in  which  edge 
averages  of  the  flux  are  approximated  with  second  order  accuracy.  It  has  to  be  noted  that, 
as  usual,  this  means  that  the  resulting  scheme  formally  only  has  a  first  order  truncation 
error.  In  the  same  way  the  first  order  schemes  used  in  Section  3  formally  only  have  zeroth 
order  truncation  error  but  it  is  well  known  now  that  they  are  convergent  (see  for  example  the 
preprints  of  Szepessy  entitled,  “Measure  valued  solutions  to  conservation  laws  with  boundary 
conditions”,  of  Champier,  Gallouet  and  Herbin  entitled  “Convergence  of  an  upstream  finite 
volume  scheme  on  a  triangular  mesh  for  nonlinear  hyperbolic  equations”,  and  of  Coquel  and 
LeFIoch  entitled  “The  finite  volume  method  on  general  triangulations  converges  to  general 
conservation  laws”). 

As  for  the  one  dimensional  case,  we  give  a  general  result  assuming  a  positive  reconstruc¬ 
tion. 


9 


We  now  consider  an  approximation  of  the  two  dimensional  Euler  equation  (3.1)  under 


the  form 


A  /  3 

up'  =  u*  -  \F\  £  hWUu"K°’,/°M  «-D 

where  we  use  the  notations  of  Section  3.  We  still  assume  that  the  approximate  Riemann 
solver  h  satisfies  the  one  dimensional  positivity  property  Assumption  1.  The  main  difference 
with  the  situation  in  Section  3,  is  that  we  now  use  second  order  approximations  U’^a ,  f/jV?  ^ 
of  Un  at  the  center  of  the  edge  e£-.  U^'a  is  such  an  approximation  in  the  triangle  K ,  is 

in  the  triangle  K(a )  neighboring  K  along  the  edge  e£-.  These  values  can  be  obtained  using 
a  slope  reconstruction,  or  an  interpolation  together  with  the  interpretation  of  functions  as 
piecewise  constant  in  subcells  as  in  Perthame  and  Qiu  [7],  or  being  evolved  as  independent 
degrees  of  freedom  in  the  discontinuous  Galerkin  finite  element  method  [2].  In  any  case,  we 
can  assume  that  these  values  are  conservative,  i.e.  that  there  exists  real  numbers  (<*>£- )i<a<3, 


such  that 


"Jr  +  «*  +  “&  =  in  “K  >  0;  E  =  \KW« 


For  instance,  when  f/£’“  is  the  value,  in  the  middle  of  the  edge  e^,  of  a  linear  function 
in  A',  the  coefficients  are  just  |A"|  times  the  barycentric  coordinates  of  the  mass  center 
of  K ,  with  respect  to  the  middle  of  the  edges.  Another  example  can  be  found  in  [7]. 

In  order  to  state  our  result,  let  us  introduce  some  notations:  Denoting  by  C £  the  middle 
point  of  the  edge  e^-,  we  define 

Y  y-'  uk^k  /4  ^ 

S  i*i 

then,  the  triangle  K  is  naturally  divided  into  three  subtriangles  (Ka)i<a<3  obtained  by 
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joining  Yk  to  the  three  vertices  of  K. 


Fig.  1  Sub-triangles  decomposition  of  a  triangle  K 

The  unit  outward  normals  of  the  sub-triangle  Ka  are  denoted  by  (i/£’“  )i<a'<3  where 
i ^'a  =  i/£-;  and  the  length  of  the  edges  (e^’0'  )i<a<<3,  where  e^Q  =  of  the  sub-triangle 
hi a ,  are  respectively  (Z^’Q  )i<a'<3  with  =  Z <j!c. 

We  can  now  state  our  main  result: 


Theorem  4:  Let  h  satisfy  the  one  dimensional  positivity  property  Assumption  1 ,  and  assume 
(4.2).  Then  the  scheme  (4.1)  satisfies  the  positivity  property  under  the  CFL  condition 

3 


T.  MKI  +  CK )  ^  a0^K 


a'=l 


Proof  of  Theorem  4:  Mimicking  the  proof  of  Theorem  3,  let  us  define  U^+1,a  by 

3 


vT'^k  =  K°K  -  Ai  £  mw6,  UT,  v°/)r/ 

0=1 


where 


Ur-°  =  Km  >  K"'e  =  Kf  lor  9^  a 


Adding  the  equalities  (4.5)  for  a  =  1,2,  3  we  obtain  (4.1)  with 

K*'  =  E 

a= 1 

Indeed,  this  is  just  a  consequence  of  (4.2)  and,  for  a  ^  /?, 

KK“'\vry/)  =  'K”  A°) 

;a,0  ;0,a 

lK  ~  lK 


(4.4) 


(4.5) 


(4.6) 


(4.7) 


(4.8) 

(4.9) 
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Now,  applying  Theorem  3  to  (4.5),  we  obtain  that  th<  U^+l'a  satisfy  the  positivity  prop¬ 
erty.  By  concave  combination  U £-+1  satisfies  it  too  and  the  theorem  is  proved. 


□ 


Notice  that  (4.2)  is  not  essential.  It  allows  us  to  simplify  greatly  the  statement  of  Theorem 
4.  If  it  is  not  satisfied,  we  require  a  two  dimensional  extension  of  the  assumption  (2.10)  for 
the  one  dimensional  case  which  could  be  quite  difficult  to  check  here.  Finally,  we  refer  to 
[7]  for  an  example  of  a  reconstruction  which  satisfies  the  positivity  property,  together  with 
(4.2).  There  the  control  volume  is  a  dual  mesh  and  the  positivity  is  proved  using  Boltzmann 
solvers.  Our  approach  here  could  be  extended  to  a  dual  mesh  and  the  full  scheme,  using  a 
Lax-Friedrichs  or  Godunov  solver,  satisfies  also  the  positivity  property. 

5  Conclusion 

We  have  considered  how  to  preserve  the  positivity  of  density  and  pressure  for  solving  com¬ 
pressible  Euler  equations  using  finite  volume  numerical  methods.  A  general  framework  is 
established  to  obtain  positivity  preserving  first  and  higher  order  schemes  for  one  and  two 
space  dimensions. 
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6  Appendix:  Positivity  of  the  Lax-Friedrichs  Scheme 

We  consider  in  this  Appendix,  the  Lax-Friedrichs  scheme  for  the  two  dimensional  situation 
in  Section  3.  Fixing  a  unit  vector  v,  the  velocity  vector  can  be  written  as  u  =  (u,tt>)  in  the 
frame  (t/,  vL).  Then,  we  define 

U  =  {p,pv,pw,E),  E  =  ^p{v2  +  w2)  +  pe,  (6.1) 

H(U)  =  F(U)  •  v  =  (pv,  pv2  +  p,  pv w,  ( E  +  p)v )  (6.2) 

The  Lax-Friedrichs  flux  is 

*(V,V)  =  i(ff(t/)+//(V)-/»(V-t/))  (6.3) 

with 

(9  =  IKN  +  «)IL  (6-4) 

Let  us  consider  f/",  f/"±1  with  positive  density  and  pressure,  and  set 

(/;+'  =Uj  -  X  [MC"h,C”)  -  (6.5) 

Then  f/"+1  also  has  positive  density  and  pressure.  An  easy  way  to  see  that  (Sanders  [8]) 
consists  of  introducing  a  splitting  of  the  equation 

Ut  +  H(U)X  =  0  (6.6) 

by 

Ut  +  (//((/)  +  0U)X  =  0  (6.7) 

Ut  +  (//((/)  -  pU)x  =  0  (6.8) 

In  other  words,  we  write  an  exact  solver  for  (5.7)-(5.8),  which  is  easy  because  there  are  no 
longer  wave  interactions  with  the  choice  of  0  in  (6.4): 

u  =  Uj-x  [«((/«)  +  dw;  -  «({/;., )  - 

0  =  Uj  -  a  [//((/■;, )  -  » (/”+1  -  H(Oj)  + 
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(6.9) 

(6.10) 


Now  (/;i+1  =  T-  indeed  the  solution  to  the  Lax-Friedrichs  scheme  (6.5),  and  since  the 
exact  solver  preserves  positivity  (notice  that  the  addition  of  ±/3U  is  unessential  for  positivity 
by  Galilean  invariance),  by  concave  combination  f/,n+1  satisfies  also  the  positivity  property. 
We  need  the  CFL  condition 

a||(M  +  c)IL<5  (6.11) 

in  order  to  avoid  the  interaction  of  waves  when  solving  (5.7)-(5.8)  by  (5.9)- (5. 1 0) . 

Another  version  of  Lax-Friedrichs  scheme  consists  of  introducing  exact  solvers  on  a  stag¬ 
gered  grid.  It  will  also  obviously  satisfy  the  positivity  property.  Both  of  these  Lax-Friedrichs 
schemes  in  addition  satisfy  the  maximum  principle  on  the  specific  entropy  (Tadmor  [10],  [5]). 

We  would  like  to  conclude  this  Appendix  by  a  remark  raised  in  [7]  on  two  dimensional 
schemes.  Usually  two  dimensional  schemes  have  to  be  written  under  the  form  (3  4)  where  the 
approximate  solver  h(U,  V ,  i/)  is  indeed  a  function  of  three  parameters  U,  V  and  v,  satisfying 
h(U,  U,  v)  =  F{U)-v.  In  [7],  the  notion  of  “genuinely  multidimensional  solver”  is  introduced, 
where  the  approximate  solver  is  indeed  under  the  special  form  h(U,V)-is  satisfying  h(U,  U)  = 
F(U).  The  Lax-Friedrichs  scheme  is  obviously  not  “genuinely  multidimensional”  because  its 
value  really  depends  on  the  frame  {u,  t'x)  used  to  write  the  four  components  system  (6.1)- 
(6.2). 

The  only  example  we  know  of  a  “genuinely  multidimensional”  solver  is  a  particular  class 
of  Boltzmann  solvers  introduced  in  [7]. 
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